d = fread("wdi.csv") %>% na.omit()
View(d)
d$SI.POV.GINI = d$SI.POV.GINI/100
d$gini
d = fread("wdi.csv") %>% na.omit()
d = d[energy_pc_gj> 50]
merge = rbind(
data.table(Gini = mod$pay_gini, Adaptation = mod$Adaptation, Type = "Energy-Hierarchy Model"),
data.table(Gini = k$Gini, Adaptation = k$Adaptation, Type = "Empirical"),
data.table(Gini = d$SI.POV.GINI, Adaptation = "Industrial", Type = "Empirical")
)
merge = rbind(
data.table(Gini = mod$pay_gini, Adaptation = mod$Adaptation, Type = "Energy-Hierarchy Model"),
data.table(Gini = k$Gini, Adaptation = k$Adaptation, Type = "Empirical"),
data.table(Gini = d$gini, Adaptation = "Industrial", Type = "Empirical")
)
merge$Adaptation = reorder(merge$Adaptation, merge$Gini, FUN = median)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_lenski.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
blank = ggplot() + theme_void()
grid.arrange( arrangeGrob(gA, gB, ncol = 2),
blank,
arrangeGrob(gC, gD, ncol = 2),
heights = c(20,1,20)
)
dev.off()
png("mod_empirical_scatter.png", width = 7, height = 6.2,  units = 'in', res = 600)
grid.arrange( arrangeGrob(gA, gB, ncol = 2),
blank,
arrangeGrob(gC, gD, ncol = 2),
heights = c(20,1,20)
)
dev.off()
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
wd = gsub("Figures", "Empirical Data/Pre Industrial", dir)
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(ineq)
dir = here()
text.size = 10
# energy-hierarchy model
##################################################
# model data
wd = gsub("Figures", "Hierarchy Model/results", dir)
setwd(wd)
mod = fread("inequality.csv")
names(mod) = c("firm_size", "power_gini", "pay_gini", "pay_exponent", "power_top_1", "pay_top_1")
# GEM firm size data
wd = gsub("Figures", "Empirical Data/GEM", dir)
setwd(wd)
energy = fread("energy_v_firm.csv")
# regression
x = log(energy$energy)
y = log(energy$firm_mean)
r = lm(y ~ x)
a = exp(coef(r)[1])
b = coef(r)[2]
e_predict = function(x){
y = (x/a)^(1/b)
return(y)
}
energy_predict = e_predict(mod$firm_size)
mod$energy = energy_predict
# Kohler with energy estimates
####################################################
# energy bounds
wd = gsub("Figures", "Empirical Data/Energy Adaptation", dir)
setwd(wd)
bounds = fread("energy_boundaries.csv")
# kohler
wd = gsub("Figures", "Empirical Data/Kohler", dir)
setwd(wd)
kohler = fread("kohler.csv")
# merge with energy bounds
kohler_energy = merge(bounds, kohler, by = "Adaptation")
# plot
kohler_plot = ggplot() +
geom_point(data = mod, aes(x = energy, y = pay_gini, color = pay_exponent), size = 0.01, alpha = 0.2) +
geom_errorbarh(data = kohler_energy, aes(xmin = `energy_5%`, xmax = `energy_95%`, y = Gini), size = 0.2, height = 0) +
geom_point(data = kohler_energy, aes(x = `energy_50%`, y = Gini), size = 0.8) +
scale_x_log10("Energy Use per Capita (GJ/year)",  breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000)) +
scale_y_continuous("Gini Index", breaks = seq(0,1, 0.1)) +
scale_color_gradientn(expression(beta), colours = rainbow(8), breaks = seq(0, 1, 0.1)) +
coord_cartesian(xlim = seq(2.5, 2000), ylim = c(0.08, 0.8)) +
ggtitle("A.  Ancient Societies") +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
legend.position = "none",
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust= 0, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
# milanovic
#####################################
wd = gsub("Figures", "Empirical Data/Pre Industrial", dir)
setwd(wd)
milanovic = fread("milanovic_energy.csv")
wd = gsub("Figures", "Empirical Data/Pre Industrial", dir)
setwd(wd)
milanovic = fread("pre_industrial_ineq_energy.csv")
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Pre Industrial/gdp_energy.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
kcal_to_gj = 0.000004186798
top_frac = 0.01
frontier_gini_func = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_subsistence_low = kcal_to_gj*2000*365
energy_low = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
top_low = mapply(top_share, energy_subsistence_low, energy_low, top_n)
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
top_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(ineq)
energy_subsistence_high = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
top_high = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(ineq)
text.size = 10
# inequality frontier GINI
####################################################################################
frontier_gini_func = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
# inequality frontier with low subsistence energy
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
frontier_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
# inequality frontier with high subsistence energy
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
frontier_high = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
frontier = rbind(
data.table(energy = energy_low, top = top_low),
data.table(energy = energy_high, top = top_high)[order(-top)],
data.table(energy = energy_low, top = top_low)[1,]
)
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
frontier_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
energy_subsistence_high = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
frontier_high = mapply(frontier_gini_func, energy_subsistence_high, energy_pc, top_frac)
energy_subsistence = 3
energy_max = 1000
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
}
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfuction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
plot(frontier_low)
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, gini = gini_vec)
return(output)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
plot(frontier_low)
plot(frontier_low, log = "x")
energy_subsistence_low = kcal_to_gj*3000*365
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max
energy_max = 3000
energy_subsistence_low = kcal_to_gj*3000*365
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*3000*365
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
energy_subsistence_low = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*4000*365
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_gini_func(energy_subsistence_high, energy_max, 0.01)
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_gini_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-gini)],
frontier_low[1,]
)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
gini_frontier
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
gini_frontier
frontier_top_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
top_share = energy_top/energy_bottom
return(top_share)
}
top_share_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, gini = top_share_vec)
return(output)
}
# calculations
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*2500*365
energy_subsistence_high = kcal_to_gj*4000*365
frontier_low = frontier_top_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_top_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-gini)],
frontier_low[1,]
)
frontier_top_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
top_share = energy_top/energy_bottom
return(top_share)
}
top_share_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, top_share = top_share_vec)
return(output)
}
# calculations
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*2500*365
energy_subsistence_high = kcal_to_gj*4000*365
# inequality frontier with low subsistence energy
frontier_low = frontier_top_func(energy_subsistence_low, energy_max, 0.01)
# inequality frontier with high subsistence energy
frontier_high = frontier_top_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-top_share)],
frontier_low[1,]
)
frontier = rbind(
frontier_low,
frontier_high[order(-top_share)],
frontier_low[1,]
)
top_frontier = ggplot() +
geom_polygon(data = frontier, aes(x = energy_pc, y = top_share), fill = "black", alpha = 0.2, col = "black") +
geom_point(data = mod, aes(x = energy, y = power_top_1), size = 0.2, col = "grey30") +
scale_x_log10("Energy Use per Capita (GJ/year)",  breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000)) +
scale_y_continuous("Gini Index", breaks = seq(0,1, 0.1)) +
scale_color_gradientn(expression(beta), colours = rainbow(8), breaks = seq(0, 1, 0.1)) +
coord_cartesian(xlim = seq(2.5, 2000), ylim = c(0, 1)) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
legend.title = element_text(hjust = 0.2),
legend.key.height = unit(1.2, "cm"),
legend.justification = c(0, 0.3),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
top_frontier
View(frontier_low)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
top_frontier
gA = ggplotGrob(gini_frontier)
# export
####################################################
setwd(dir)
gA = ggplotGrob(gini_frontier)
gB = ggplotGrob(top_frontier)
png("inequality_frontier.png", width = 7.5, height = 4,  units = 'in', res = 600)
grid.arrange(cbind(gA, gB, size = "first"))
dev.off()
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_result.R')
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(gtable)
text.size = 10
# model
#################################################################
dir = here()
wd = gsub("Figures", "Hierarchy Model/results", dir)
setwd(wd)
mod = fread("inequality_energy.csv")
plot(mod$alpha, mod$power_gini)
plot(mod$alpha, mod$power_gini, rev = "x")
plot(mod$alpha, mod$power_gini)
plot(mod$alpha, mod$power_gini)
plot(mod$energy, mod$alpha)
plot(mod$energy, mod$alpha, log = "x")
plot(mod$energy, mod$alpha, log = "xy")
plot(mod$energy, mod$alpha-2, log = "xy")
plot(mod$energy, mod$alpha, log = "xy")
plot(mod$energy, mod$alpha, log = "x")
plot(mod$alpha, mod$energy, log = "x")
plot(mod$alpha, mod$energy, log = "y")
plot(mod$alpha, mod$pay_gini)
